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Abstract: The paper address the reahzation and identification problem for a subclass of 
piecewise-affine hybrid systems. The paper provides necessary and sufficient conditions for 
existence of a realization, a characterization of minimality, and an identification algorithm for 
this subclass of hybrid systems. The considered system class and the identification problem are 
motivated by applications in systems biology. 
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1. INTRODUCTION 



In this paper we address the realization and identification 
problem for a subclass of piecewise-affine hybrid systems. 

Contribution of the paper We define the class of 
piecewise-linear systems (abbreviated by PWL). PWLs 
are a subclass of piecewise-affine hybrid systems. The 
continuous dynamics of PWL is determined by a finite 
collection of affine subsystems. However, in contrast to 
traditional piecewise-affine systems, we allow any change 
of the continuous state during a discrete-state transition, 
as long as the new state belongs to the set of designated 
initial states of the aSine subsystem associated with the 
new discrete state. In addition, we do not impose any spe- 
cific mechanism for triggering discrete-state transitions. 

We formulate the realization problem for this system class 
and partially solve it by providing n(x;cssary and sufficient 
conditions for existence of a realization. We also present 
conditions for minimality. We show that the outputs 
of any PWL can also be described by a switched AR 
model. The main conclusion is that any system can be 
transformed to a minimal system with one discrete state 
while preserving input-output behavior. This means that 
without further restrictions, the identification problem for 
such systems is not necessarily interesting. We discuss a 
number of restrictions on the system structure which avoid 
this problem. 

Note that the conclusion above is not valid for other classes 
of hybrid systems. For hybrid systems from Petreczky 
and van Sc;huppen (2010), there examples of input-oiitput 
maps which provenly cannot be realized by systems with 
one discrete state. 

In addition, we present an identification algorithm for sys- 
tems with full observations. This algorithm is illustrated 
by example of physical and biological relevance. 

Motivation The motivation for studying realization the- 
ory for PWLs is that it provides the theoretical founda- 



tions for systems identification. The motivation for inves- 
tigating identification of PWLs is the following. 

• First, several systems of interest can be modeled by 
PWLs and PWLs tend to be convenient for analysis. 
However the parameters of PWL models are often 
not directly available and hence they have to be 
estimated from measurements. 

• The problem of estimating the network dynamics of 
complex biological systems can be reduced to the 
identification problem for PWLs. 

Identification of PWLs and biological netw^orks 
Below we will elaborate on the relationship between iden- 
tification of PWLs and the estimation of dynamic inter- 
actions of complex dynamical systems. 

Numerous vital processes in nature involve complex signal- 
ing networks, varying from gene-protein interaction net- 
works Westra et al. (2007) and complex communication 
in microbes Gera and Srivastava (2006) to the synchro- 
nization in the heart of higher animals Heijman et al. 
(2009). In most cases, these signaling networks do not 
have a fixed topology, but vary their structure according 
to their internal states and certain external conditions 
Adami et al. (2000). This flexibility allows the organism to 
exhibit complex behavior and flexible responses to various 
environmental conditions. 

The underlying biological problem is to describe these flex- 
ible dynamical networks, based on (partial) observations 
of their internal states. Examples of potential observations 
are time series of gene expressions, protein densities, oxy- 
gen stress and sugar concentrations. 

In almost all cases, there is no explicit mathematical 
model available that even remotely describes the network 
dynamics. This is caused by the sheer size of the problem; 
the many thousands of genes, RNAs, proteins and other 
types of molecules involved, the complexity and idiosyn- 
crasy of their interactions, and the large amount of noise 
in biological systems. Basically, many agents involved in 
the network are unknown let alone their mechanism of 
interaction. 



Mathematically, the processes above can be viewed as 
nonlinear dynamical systems with unknown parameters. 
Such systems often have several equilibria, and hence their 
behavior can be seen as a mixture of certain afRne systems, 
where each afRne system is obtained by linearizing the 
original system around one of the equilibria. That is, such 
systems can be approximated by a PWL. In addition 
to their simplicity, the advantage of PWL approxima- 
tions is that they neatly capture the interaction among 
various state components around various equilibria. More 
precisely, interaction between the ith and jth state compo- 
nents can be viewed as the property that the (z, j)th entry 
of the Jacobian of the system around equilibrium is non- 
zero. Since we are often interested exactly in interactions 
rather than the detailed dynamics, it makes sense to recast 
the problem of identifying such interactions into the prob- 
lem of identifying PWLs approximations. While the paper 
does not directly address identification of interactions, we 
believe that the presented results represent a step towards 
the solution of that problem. 

Related work Identification of piecewise-linear (-afiine) 

hybrid systems has been subject of intensive research, 
Vidal et al. (2002); Ma and 'Vidal (2005); Bako et al. 
(2009b); Paolctti ct al. (2007); Vidal (2008); Ferrari- 
Trecate et al. (2003); Juloski et al. (2004); Bako et al. 
(2009a); Roll et al. (2004); Fox (2009). Realization theory 
for hybrid systems was investigated in Weiland ct al. 
(2006); Paoletti et al. (2010); Petreczky (2006); Petreczky 
and van Schuppen (2010); Grossman and Larson (1995). 
The application of hybrid systems to modelling complex 
dynamical systems, originating from biology, and their 
identification has been a subject of intensive research 
Casey et al. (2005); de Jong (2002); Porreca et al. (2010); 
Koutroumpas et al. (2007); Cinquemani et al. (2008). 

The results of this paper are new, to the best of our 
knowledge. What distinguishes the contribution of this 
paper from the existing work is (a) the system class 
considered, in particular, the emphasis on continuous- 
time systems, (b) the emphasis on realization theory and 
algorithms, (c) the details of the identification algorithm. 

Outline of the paper §2 presents the formal definition of 
the system class of interest. §3 discusses the relationship 
between nonlinear systems with complex dynamics and 
picccwise-affinc hybrid systems. §4 presents the results 
on realization theory, and §5 presents the identification 
algorithm and the results of the numerical experiments. 

Notation We use the standard notation. We denote by 
T = [0,-|-oo) the time-axis. We denote by 7„ the n x n 
identity matrix. We denote the set of natural numbers 
including zero by N. 

2. PIECEWISE-LINEAR SYSTEMS (PWL) 

The aim of the section is to define the class of piecewise- 
hnear systems formally. 

Definition 1. (PWL). A piecewise-linear system (abbrevi- 
ated as PWL) is a dynamical system determined by 



Here Q = {1, . . . , D} is the finite set of discrete modes, MP 
is the output space, M"' is the state-space of the system 
in mode q € Q, For each q,G Q Ag G M"? 



hit) 



aq G 



, and Cq e are the parameters of the affine 
system in mode q G Q, Xq^ C M" - is the set of initial 

states of the affine system in mode q € Q. The state space 



ofSis7^E=U,eQW 



We call S linear, if c„ = 



y{t) = Cg(t)X{t) + Cq(t) 



(1) 



and ttq = for all g € Q, otherwise S is called affine. 
We will use the following short-hand notation 

S = Q, {{ng, Ag, Ug, ^ Q ) l ^ ^ ^1). 

Informally, the evohition of E takes place as follows. As 
long as the value of the discrete state q{t) at time t 
does not change, the continuous state and the continuous 
output an time t change according to the affine system 
i{t) = Aqi^t)x{t) + aq(t) and y{t) = Cq(t)x{t) + Cq(ty The 
discrete state can change at any time, however, we do not 
allow consecutive changes of discrete states immediately 
one after the other. If the discrete state changes to q{t'^) 
at time t, then the new continuous state should satisfy 
x(t+) e A'g(t+),o: i-t!- it should belong to the set of the 
designated initial states of the new discrete state qif^)- 
Note that we do not specify the mechanism which triggers 
the change of discrete state. We also do not specify 
the initial discrete state, it is chosen by an unspecified 
mechanism. Hence, the description above allows for several 
state- and output-trajectories. 

Note that the definition of PWLs presented above differs 
from the one used in the literature. In the standard defi- 
nition it is assumed that after a discrete-state transition, 
the new continuous state depends on the previous one. 
In contrast, here we only require that the new continuous 
state belongs to the designated set of initial states of the 
affine system associated with the new discrete state. This 
implies that the behavior of the system in discrete state 
q depends only on the affine system associated with that 
discrete state and it does not depend on the discrete modes 
which were visited in the past. Another consequence of the 
definition is that the model is not predictive. That is, the 
knowledge of the system parameters and the inputs does 
not uniquely determine the state- and output-trajectory 
of the system. 

In order to define the evolution of PWLs formally, we 
need to introduce the following notation and terminology. 

Definition 2. (Collins (2005)). A time event sequence is a 
strictly monotone sequence {tn)n=o ^^ch that n* G N U 
{+00}, to = and for all < n < n*, < t„ < t„+i. If 
n* = -l-oo then we require that sup{i„ | n G N} = -|-oo. If 
n* < -foo, then by convention tn*+i — +00. 

The role of time event sequences is to formalize the time 
instances at which discrete events occur. The restrictions 
formulated in the definition imply that the set of switching 
times does not have an accumulation point, i.e. no Zeno- 
behavior can take place. 

Definition 3. (State-trajectory). A state-trajectory of E 
is a map ^ : T — >■ 'Hy. such that there exists a time 
event sequence (ti)"^o ''^^'^ ^ sequence of discrete modes 
(ft G <9)"=Io such that for all < i < n*, i G N, it holds 
that for all s G [fi,ti+i), ^(s) = {qi,x{s — tt)) and 

x{t) = Ag.x{t) + Gq^ and a;(0) G Xq.^. 



The time event sequence {tn)n=Q called the sequence of 
switching times of the state-trajectory ^. 

We denote by BS{Ti) the set of all state-trajectories of S. 

Definition 4- (Output-trajectory). An output-trajectory of 

S is a map j/ : T — > such that the following holds. There 
exists a state-trajectory ^ of S, such that y{t) = vs{x{t)) 
for all t gT. Here, vs is the readout-map of S, defined as 

: 3 {q,x) CqXeW. 

We denote by B{T,) the set of all output-trajectories ofT,. 

In the sequel, unless stated otherwise, / denotes a function 
f : T ^ MP. The definition above implies that the external 

behavior of a PWL is exactly a function of this type. 

Definition 5. (Realization). The function / is said to be 
realized by PWL E, if i/ is an output-trajectory of S, i.e. 
if / G -B(S). In this case S is called a realization of /. 

Definition 6. The dimension of a PWL I] is a tuple 
(|(3|,n), where n = J2qeQ^'i ^ linear, and n = 
Y^q^qiriq + 1), if S is affine. 

That is, the first element of dim S is the number of discrete 

states, the second clement is the number of continuous 
state components. The additional dimension in the case 
of afHne PWLs stems from the need to store the vectors 
ttq, Cg, q e Q. We use the following partial order relation 
on N X N. We say that (p, g) e N is smaller than or equal 
(r, s) € N, denoted by {p, q) < (r, s), ii p < r and q < s. 
Note that the order relation < in N x N is indeed a partial 
order, it is not possible to compare all elements of N x N. 

Definition 7. (Minimality). A PWL S is said to be a 
minimal realization of /, if for any PWL S which is a 
realization of /, dimS < dimE. 

The definition above implies that if S is minimal, then it 
is possible to compare its dimension to the dimension of 
any other realization of /. Since we work with a partial 
order on dimensions,the existence of a minimal system is 
not at all obvious. 

Problem 1. (Realization problem for PWLs). Find condi- 
tions for existence of a PWL realization of /. Characterize 
minimal PWLs realizations of /. Find algorithms for 
computing a PWLs realization from finite data. 

As it was indicated before, our motivation for studying 
the realization problem for PWLs is to lay the theoretical 
foundations for identification of PWLs. We present below 
the formulation of the identification problem. 

Problem 2. (Identification problem). Assume that the value 
of / and its derivatives up to order r, r > is measured at 
time instances t\ < . . . < tk- Based on the (possibly noisy) 
data {/^'^(ti)}j=i...,fe,i=i,...,ri find a PWL realization of /. 

3. APPROXIMATION BY PWLS 

Our main motivation for studying identification of PWLs 
is the following. The problem of identifying interaction 
networks of complex dynamical systems can be reduced 
to the identification problem of PWLs. Indeed, consider 
a partially observed non-linear system of the form 

x{t) = f{x{t)) and y{t) = h{x{t)) (2) 



where / : and h : ^ W are sufficiently 

smooth functions. 

We are interested in the dynamics of the network structure 
of (2). By the network structure of (2) at time t we mean 
the directed graph with nodes numbered by 1, . . . , and 
where an edge goes from node i to node j if the ith 
component of x{t) influences the )th component of x{t). 
As the values of x{t) changes, so does f{x{t)). Hence the 
network structure described above depends on time. 

The; basic question of interest is how to derive the dynam- 
ics of the interaction network based on the observed output 
y. One obvious approach is to estimate / and h from y, but 
this is feasible only if special restrictions are put on / and 
h. Notice, however, that if we look at local linearizations 
of / around an equilibrium and the interactions are strong 
enough, then Xi{t) influences Xj{t), if the {i,j)th entry of 
the Jacobian of / is non-zero. 

Prompted by the considerations above we propose the 

following approach to modelling and identification of the 
network structure of systems of the form (2). We assume 
that (2) has finitely many equilibrium points. We assume 
that under influence of noise the system jumps between 
neighborhoods of these equilibrium points. We assume 
that the transition from the neighborhood of one equi- 
librium point to the neighborhood of another equilibrium 
point takes place very fast. Hence, the behavior of the 
system during the transition can be ignored. We associate 
a PWL S of the form (1) with (2) as follows. Let D be 
the mimbcr of equilibrium points. For each q E Q, let 
be the corresponding equilibrium point and define 

Ag = Df{eq) and a, = -AqCg 
Cq = Dh{eq) and Cq = h{eq) 

Here Df and Dh denote the Jacobian of / and h respec- 
tively. Let Xq Q be a suitably chosen neighborhood of Sq. 

It then follows that if the state x{t) of (2) is close enough 
to Cg, we can approximate (2) as follows 

x{t) = f{x{t))^Aqx{t) + aq 

y{t) = h{x{t)) « CqX{t) + Cq 

As it was noted above, the relevant information regarding 
the interactions of state variables is already included in 
Aq. Hence, for the purposes of network reconstructions, 
identification of the PWLs approximation is sufficient. 

4. REALIZATION THEORY OF PWLS 

Below we present some basic results on realization theory 

of PWLs. First we present a characterization of mini- 
mality, after that we present conditions for existence of 
a realization and a realization algorithm. Throughout the 
section, / denotes a function f -.T^MP. 

4.1 Minimality 

In order to present our first result, we introduce the notion 

of a linear system with state-jumps. 

Definition 8. A linear system with state jumps, abbre- 
viated as (LSSJ) is a linear PWL E of the form (1) 
with one discrete state, i.e. D = 1. We will identify the 
LSSJ S with the collection of data (n, C, A, A'o), where 
n = ni, C = C\, A = Ai and X\fi = Xq. 



The reason we call a PWL with one discrete state a 
linear system with state jumps is that the system behaves 
as a linear system, with the exception that its state 
occasionally jumps back to one of the initial states. For 
LSSJs, we can easily define the concepts of observability 
and span-reachability. 

Definition 9. A LSSJ S = {n,C, A, Xq) is called observ- 
able, if (C, A) is an observable pair, and S is called span- 
reachable, if SpanjA'^xo \ k = 0, . . . ,n — I.xq E Xq} = n. 
Remark 1. Note that Z = Spanjyl'^xo | k — 0, ...,n — 
l,xo G Ao} is an A invariant subspace containing Xq, 
hence by restricting A and C to Z we can transform S to 
a span-reachable LSSJ while preserving output trajec- 
tories, i.e. B{'E,r) = B{Yi). By applying linear observability 
reduction, we can transform the LSSJ to a LSSJ 
S„ which is span-reachable, observable and has the same 
output trajectories as S, i.e. = 
Theorem 1. (Minimality). Assume that / admits a real- 
ization by a PWL. Then / admits a minimal PWL real- 
ization E such that E is a span-reachable and observable 
LSSJ. A LSSJ realization of / is minimal if and only if 
it is span-reachable and observable. 

The theorem above says that in general, the external 
behavior of any PWLs can be represented by a linear 
system with several initial states. That is, without further 
restrictions, the realization and identification problems for 
PWLs are equivalent to that of LSSJs. 

The proof of Theorem 1 relies on the following transfor- 
mations. 

Definition 10. (PWL to linear PWL). Define the linear 
PWLs i(E) associated with an affine E as follows. 



L(E) = (p,Q,{«,A^ 



where = -|- 1, 

'A, 



I ) ) ^q 
L - 



= 0, ai' = and 
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{(a;^,i)^eM"''+Mxe ^-,,0}. 

Proposition 1. The output trajectories of E and £(E) 
coincide, i.e. S(E) = B(I-(E)), and dimE = dimL(E). 



Proof. [Sketch] A state trajectory ^ of i(E) is always a 
map of the form ^(t) = {q{t), (a;(t)^, 1)^), t e T such that 
^{t) = {q{t),x{t)), t e T is a state-trajectory of E. Since 
^^L(E)(('7, (a;^, 1)^)) = CqX+Cq = ((s, 2;)) , the Statement 
of the proposition follows. 

Definition 11. (Linear PWL to LSSJ). Let E be a linear 
PWL of the form (1). Define the LSSJ LS'(E) associated 
with E as follows. 

LSiE) = {n,C,A,Xo), 



where n = J2q(^Q '' 
[Ai • 
Ao ■ 



and 






and C=[Ci ■■■ Cd] 



••• Ad J 

'^0 = -^9,0 

Xg,o = {{ 0^_^ ,x,0,...,Of \ x€Xg,o},yq^Q 

q—l — times 



Proposition 2. The output trajectories of E and LS'(E) are 
the same, i.e. -B(E) = B{LS{'E,j), moreover, dimL5'(E) = 
(l,n) where (|(3|,n) =dimE. 

Proof. [Sketch] For each q G Q, let Xg be the subset of 
of the form z = ( 0,0, -^..,0 , a;"^, 0, . . . , 0)^, x e M"' . For 

q—l— times 

each z G Xqof the above form, let 11^(2;) = x G M"« . It then 

follows that ^(t) = (1, x{t)) is a state-trajectory of _L5(E) if 
and only if ^(t) = {q{t),Ilq(^t){x{t))) is a state-trajectory of 
E, where q{t) = q for some q G Q if and only if x{t) G Xq. 

Since uls(s)((1, a^)) = Cqllq{x) = vj:{{q,nq{x)) for all 
X Cz Xq, the statement of the proposition follows. 

Using the two transformations defined above, and Remark 
1, we can present the proof of Theorem 1. 

Proof. [Proof of Theorem 1] Assume that E is a PWL 
realization of /. If E is affine, then replace E with the 
associated linear PWL i(E), which is also a realization of 
/. Hence, we can assume that E is already a linear PWL. 
Construct then the LSSJ ^^(E) associated with E and 
apply Remark 1 to LS{T,) to obtain a span-reachable and 
observable LSSJ E„ such that B(E„) = B(LS'(E)) = 
S(E). It then follows that dimE„ < dimLS'(E) < dimE. 

Hence, since E was arbitrary, it is enough to look for 
minimal realizations among span-reachable and observable 
LSSJs. Notice that for any two LSSJs Ei and E2, either 
dim El < dimE2 or dimE2 < dimEi, hence, among all 
the possible LSSJs realizations of /, there must exist a 
minimal one. 

4.2 Existence of a realization 

The conditions for existence of a realization will be formu- 
lated using the rank of the Hankel-matrix of /. In order to 

define the Hankel-matrix of /, we have to define the notion 
of Markov-parameters. To that end, we need the notion of 
piecewise- analytic functions. 

Definition 12. (Piecewise- analytic). The map / is called 
piecewise- analytic, if there exist a finite or infinite number 
of time instances, ti E T, ti < ti+i, i < Nf, i G N for 
some A^/ e N U {-|-cxd}, such that the following holds. 
If Nf < +00, let tN^+i = +00. Then we require that 

to = and [jfjo[ti,ti+i) = T and for each i e N, i < Nf, 
f is analytic on [ii,ti+i), but / is not analytic on any 



the 



neighborhood of ti in M. We call the points {<i}^o 
points of non-analyticity. We define the set 

If = {zeN\z<Nf} 

of indices of points of non-analyticity. 



The intuition behind the definition is as follows. If / has 
a realization by a PWL E, then the only points where 
/ is not analytic are the points where the corresponding 
state-trajectory of E switches from one discrete mode to 
another. In fact, one can show that there always exists a 
state-trajectory of E which yields / as output trajectory 
and which switches only at time instances at which / is not 
analytic. Hence, if / has a realization by a PWL, then / 
is piecewise-analytic and the points of non-analyticity tell 
us the switching times of a PWL realization of /. 

In the sequel, / is assumed to be piecewise-analytic. 



Definition 13. (Markov-parameters). Assume that / is 

piecewise-analytic and let {ti}^!^ be the points of non- 
analyticity of /. For each iGlf, define the ith Markov- 
parameter of / as a sequence : N — )• ffi*' 

VfceN:Mf(fc) = ^/(i, + s)U=o. 

It is easy to see that the collection of Markov-parameters 
{M{}^^q determines the map / uniquely. We use the 

Markov-parameters to define the Hankcl-matrix of /. 

Definition 14- (Hankel- matrix). We define the Hankel- 
matrix Hf of / as the infinite matrix, rows of which 
are indexed by N x {1, . . . and columns of which are 
indexed by N x The entry of Hf indexed by row index 
{i, r) and by column index (j, I) equals 

[Hf]{i,rUo,l) = {M.l{i+3))r 

where (M/ {i + j))r denotes the rth entry of Markov- 
parameter yif (*+.?)• The rank of Hj, denoted by rank H f\, 
is the dimension of the linear space spanned by the 
columns oi Hf. 

Theorem 2. (Existence of a P WL realization) . The map 
/ can be realized by a PWL if and only if it is piecewise- 
analytic and rank Hf < +oo. Moreover, a minimal LSSJ 
realization of / of dimension (l,ranki?/) can be con- 
structed from Hf . 

Proof. [Sketch] only if Assume that / can be realized 
by a PWL. Then it can be realized by a LSSJ S = 

{n,C,A,Xo)- Then M/(fc) = C^^x, Vi G //,fc G N for 
some Xi G Xq, i.e. {-^/(A:)}^q are the Markov-parameters 
of the linear system (C, A) from some initial condition. It 
then follows from linear systems theory that rank Hf < 

n < +0O. 

if Assume that n = rank Hf and fix a basis in the 
linear span of the columns of Hf. Define C as the matrix 
in this basis of the linear map which maps a column to 
the vector formed by the rows of that column indexed 
by (0, 1), ... , (0,p), in this order. Define the matrix A as 
the matrix in this basis of the linear map which maps 
the column indexed by (j, /) to the column indexed by 
{j + 1,1). Finally, let Xq be the set of coordinate vectors 
of the columns of Hf indexed by (0,1), I G If. Then 
'Sf = {n,C,A,Xf)) is a LSSJ realization of /. If S is 
another LSSJ realization of /, then from the only if part 
it follows that rank Hf = dimS/ < dimS, i.e. S/ is a 
minimal realization of /. 

In Algorithm 1 we present a Kalman-Ho-like realization 
algorithm for LSSJs and hence PWLs. To this end, we 

need the following definition. 

Definition 15. (Finite Hankel sub-matrix). Fix integers 
i?, L, M > 0. Define the set 

If = {iGlf\i< R}. 

Define the submatrix Hf^L^M,R of Hf as the matrix which 
is formed by the intersection of the rows of Hf indexed 
by the elements oi II = {0, ■•-,-£'} x {1, ... ,p} and 
the columns of Hf indexed by the elements of Jm,r = 
{0,...,M}x/f. 



Algorithm 1 

Inputs: Hankel-matrix Hf^L,M+i,R- 
Output: LSSJ Sl m h- 

1: Compute the decomposition Hj^l m+i.r = OR such 
that O G M^^><" and R G M"^'^"+i.« and rank R = 
rank O = n. 

2: Define R G M."'^-^^'' as the matrix formed by the 
columns of R indexed by elements of Jm,r- Let R G 
^nxJM,R ^ijg matrix such that the column of R indexed 
by (j, I) equals the column of R indexed by {j + 1, 1). 

3: Define ^l,m,r = {n,C,A,XQ) as follows. 

• The ith row of C equals the row of O indexed by 
(i,0), i = l,...,p. 

• The matrix A is the solution of the equation 

n = An 

• The set Xq consists of the columns of R indexed 
by indices of the form (0, 1), I G if. 



Theorem 3. (Correctness of Algorithm 1). If rank Hf < 

n, and |{M/ | i G //}| < -|-oo, then for some R > 0, the 
LSSJ returned by Algorithm 1 is a minimal realization 

of /. The conditions above hold if / has a realization by 
a LSSJ of dimension at most n and with a finite set of 
initial states. 

4.3 Ill-posedness of the realization problem 

The results on realization theory of PWLs indicate that 
the realization and identification problems for PWLs are 

in general ill-posed in the following sense. Realizability by 
a PWLs is equivalent to realizability by a minimal PWLs 
with one single discrete state. Hence, thcire seems to be no 
intrinsic way to choose discrete states based on the data. 
Since the original intention was use the combination of 
the dynamics in various discrete modes to explain complex 
dynamics, the above conclusion is an unpleasant one. 

Note that the main problem is that one can increase the 
dimension of the continuous state-space to encode discrete 
states. One way to remedy this is to place restriction on 
the dimension. To this end, we introduce the notion of 
K — N realization. 

Definition 16. {K - N). A PWL S is said to he &K-N 
PWL if S is of the form (1) and \Q\ < K and for all 9 G Q, 

Uq < N . The map / is said to admit a K — N realization, 
if there exists a K — N PWL which is a realization of /. 

The motivation for using K — N realizations is that by 
choosing an appropriate N, we can make sure that the 

growth in the number of continuous states cannot be used 
to replace discrete states. Below we present conditions for 
existence of a K — N realization of /. To this end, we have 
to introduce the following concept. 

Definition 17. (Hankel-matrices & partition). Consider the 
partitioning C = (Cg)^i of the set of If. For each 
q — 1, . . . , if, let Hj be the sub-matrix of Hankel-matrix 
Hf which is formed by the columns oi Hf indexed by 
indices of the form {j, /), j G N and I G Cq. 

Theorem 4- (Existence oi K — N realizations). The map 
/ has a realization by a linear K — N PWL if and only if 
there exists a partitioning C = (0,)^^ of I with D < K 



such that for alH = 1, . . . , D, rank Hj < TV. If the latter 
condition holds, then S can be computed from Hf with 
Q = {1, . . . , D} and with = rank Hj f-, for all q G Q. 

Theorem 4 above says that / can be realized hy a, K — N 
realization if the columns of the Hankcl-matrix Hf can 
be divided in K clusters such that for each fixed I G /, 
all columns indexed by end up in the same cluster. 
Moreover, the linear span of the elements of each cluster 
should have dimension at most N. 

In order to prove Theorem 4, we need the following 
transformation of a LSSJ to a. K — N realization. 

Definition 18. (LSSJ to PWL). Consider a partition C = 
{^i)iLi of If for some K > Q. Consider an observable 
LSSJ S; = (n, C, A.Ao) such that S; is a realization of 
/. Then there exists a unique state trajectory ^ of S;, 
such that the switching times of ^ coincide with the points 

of non-analyticity (ii)i=fo of ./ and = I{t) for 

all t e T. Define the PWL realization PW{Y.i,'C) of / 
associated with S; and C as fohows. The PWL PW{T,i, C) 
is of the form (1), such that Q = {1, . . . , if}, and 

• For all q & Q, let Xq^ be the set of states x{ti) such 
that for some I € C^, ^{ti) — {l,x{ti)). Define 

Xq = Span{A''x | fc = 0, . . . , n - 1, a; e Xq^o} 

Set Uq = dimXq and choose a basis of Xq. 

• For each q £ Q, define Aq as the matrix in the basis 
of Xq of the linear map Xq B x i-^ Ax G Xq obtained 
by restricting A to Xg. 

• Let Xqfi be the set of vector representations in the 
basis of Xq of the elements of Xq^. 

• Let Cq be the matrix representation in the basis of 
Xq of the linear map Xq B x t-^ Cx G W obtained by 
restricting C to Xq. 

Proposition 3. If S; is a LSSJ realization of /, then the 
PWL PW{Y.i, C) is a rcahzation of /. 

Proof. [Sketch of the proof of Theorem 4] only if As- 
sume that S is a PWL realization of / of the form (1), such 
that |Q| < and Uq < N for all g € Q. Then there exists a 
state trajectory ^ of E such that / = f s and the switch- 
ing times of ^ contain the points of non-analyticity {ti)^J.Q 
of /. Define Cq = {i G If \ ^{ti) = {q,Xi) for some Xi € 
Xqfi}. It then follows that C = (Cq)^i is a partitioning of 

//. Moreover, for any i e C,, M/ is the Markov-parameter 

of the linear system [Cq,Aq) from some initial condition. 
Hence, by using linear systems theory we can show that 
rank H'f q <nq < N . 

if If the conditions of the theorem hold, then rank Hf = 
n < +00. Construct the LSSJ realization Ej of /, as in 
the proof of Theorem 2. Apply Definition 18 to obtain a 
linear PWL E PVF(E/, C). Then E has D < i\: discrete 
states. Notice that the space Xq from Definition 4 is then 
isomorphic to the column space of H'j q = 1,. . . ,D. 
Hence = dim Xq = rank Hj q < N. 

Similarly to the general case, in Algorithm 2 we state an 
algorithm for computing a, K — N realization of /. 

Theorem 5. (Correctness of Algorithm 2). If / has a, K — 
N realization of dimension at most (i^, min{L, M}), and 



Algorithm 2 

Inputs: Hankcl-matrix H f ij^i^i u- 
Output: PWL E 

1: Apply Algorithm 1 to Hf^L,M+i,R and denote the 
result by T.l,m,r- 

2: Compute a partitioning C = (Ci)-^i of If such that 
D < K and for alH = 1, . . . , rank Hf c l m k — ^ ■ 
Here H^f c,l,m r ^^'^ sub-matrix of Hf^L,M,R formed 
by the columns (j, I) where j < M and I £ Cq. 

3: Return the PWL E = PW(El,m,k, C). 

and with a finite set of initial states, then for some R > 0, 
Algorithm 1 returns a minimal realization of /. 

4-4 Equivalence of PWL and switched AR models 

Below we show that / is realizable by a PWL if and only 
if / satisfies a switched AR model. 

Definition 19. (SARS models). A switched AR system 
(abbreviated as SARS) is a tuple 

I = [p, n,Q, {Ag,^ \ qeQ,i = l,...,n}) (3) 

where Q = {1, . . . , D}, D > 0, nq > and Aq^i e Rp^'p 
for alH = 1, . . . , rig. The function / is said to satisfy the 

SARS model T, if the following holds. Let (ti)^^o 
points of non-analiticity of /. For any i G If, and for any 
t G [ti,ti-^-i), denote by f^''\t) the kth order right-hand 
derivative j^f{t + .s)|s=o- Then we require that for any 
i G If there exist q{i) G Q, such that Vt G [t,,ij+i), 

/(")(t) = Xi^W,fe/^"-'=nt)- (4) 
fe=i 

The set Q is called the set of discrete modes of I. 

Notice that if D = 1 , then the SARS I corresponds to an 

AR system, and / satisfies I if it satisfies an AR equation. 
Note, however, that in contrast to trajectories of an AR 
systems, / is not smooth. 

Theorem 6. The fimction / has a realization by a PWL 
with D discrete states if and only if there exists a SARS 
I with D discrete modes such that / satisfies I. 

Proof. [Sketch] only if Let E be a realization of /. 
Without loss of generality we can assume that E is linear, 
it is of the form (1) and Uq = n for all q G Q. Then for each 
i G If, there exists q{i) G Q such that on [ti,ti^i), f is an 
output trajectory of the linear system (Cg(j), Ag(j)). From 
classical theory we then know that there exists matrices 
. . . , A5(i),„ such that (4) holds. 

if Assume / satisfies a SARS X of the form (4) . Define 
the linear PWL Ex of the form (1) such that Ex is a 
realization of / as follows. The sets of discrete modes of 
Ex and I coincide, Uq = n, Xq^ = M"" and the matrices Cq 
and Aq, q G Q are the matrices of the linear system which 
corresponds to the AR ?/(")(*) = Ei=i Here 
we consider the linear system, state vector of which is the 
regressor x{t) = ((y(""^))'^(i), • • .,y^{t)Y'. 
Theorem 1 and Theorem 6 yield the following. 
Corollary 1. The function / has a realization by a PWL 
if and only if it satisfies an AR system, i.e. for some 
Ax,...,AnGW"'P, n > 0, 



i=l 

5. IDENTIFICATION ALGORITHM FOR PWLS 

Below we present an algorithm which computes a real- 
ization with full observations based on measurements at 
finitely many time instances. More precisely, the algorithm 
solves the following problem. 

Problem 3. Fix integers D > and n > 0. Consider 
a piecewise-analytic function / : T — >■ M" and a finite 

sequence of time instances ti < ... < tM and assume 
that {,f{ti),f{ti)}f^^ are known, i.e. we know that value 
of / and its derivative at time instances ti, . . . ,tk- Find a 
PWL S of the form (1) with full-observations such that E 
realizes /, and Cg = /„, Cq = 0, and Uq = n for all q € Q. 

Motivation of the identification problem The moti- 
vation for considering the identification problem with full 

observations is the following. 

(1) Notice that the identification problem with full obser- 
vations and the identification problem for SARSs are 

equivalent. Indeed, a PWL with full observation can 
be considered as a SARSs with n = 2. Conversely, 
for a SARS I, the associated PWL Ex from the 
proof of Theorem 6 can be viewed as a PWL with 
full observations, if the high-order derivatives of / can 
be measured. Hence, by Theorem 6, the identification 
problem for PWLs is equivalent to the identification 
problem for PWL with full observations, if we assume 
that high-order derivatives of the output can be mea- 
sured too. 

(2) The problem of identification with full observations is 
still a non-trivial problem. Even in this case we have 
problems with identifiability, see Example 1. 

(3) The solution of Problem 3 enables us to prove ex- 
perimentally the feasibility of approximating complex 
systems by PWLs. This is done by applying the 
identification algorithm which solves Problem 3 to 
simulated trajectories of several well-known complex 
systems. 

Example 1. Consider the LSSJs Si = (2, Ci, Ai, A'q^) and 
S2 = (2, C2, A2, A^Q^) where Ci = C2 = h is the identity 

= {(1,0)-^} and the remaining 



matrix and Xq = 



parameters are as follows; Ai 




1 



and ^9 



2 
3 



It then follows that f{t) = (1,0)-^, t £ T can be realized 
both by Si and S2, but clearly Si 7^ S2. In other words, 
realizations of f with full observation are not identifiable. 
The reason for this phenomenon is that the state of both 
Si and S2 live in the subspace (x,y), y = 0, hence the 
difference in and A2 is not visible. 

Remark 2. (Related work on SARSs). There is a wealth 
of results on identification of SARS, mostly addressing 
the discrete-time SISO case. Many of these results, in 
particular, the algebraic approach Ma and Vidal (2005); 
Vidal (2008); Bako et al. (2009b) could probably be 
adapted to the system class of this paper. Investigating 
such adaptations remains a topic of future research. Note 
that it is not at all clear that the convergence results for 
the discrete-time case have meaningful counterparts in the 
continuous-time case. 



In Algorithm 3 we present a solution for solving Prob- 
lem 3. Algorithm 3 is an iteration consisting of the fol- 
lowing steps. First, the algorithm starts with a random 
initialization of {Aq,aq}q^Q. Subsequently, in each itera- 
tion step, the current estimates ^9, flq, {Wg,i}qgQ,i=i,...,M 
arc updated by applying first the step 3 for estimating 
the updated weights with the fixed system parameters 
{Aq, aq}q^Q, and then applying step 4 with the previously 
updated weights {w^g,i}9eQ,j=i,...,M for updating the esti- 
mates of the system parameters {Aq,aq}q^Q. The inter- 
pretation of the weights is as follows: Wq i is one if in tj 
the discrete state q € Q is active, and it is zero otherwise. 
The algorithm terminates when an absolute criterion E 
falls below a pre-specified threshold e. The iteration fails 
if after a pre-defined maximum number of iterations T^ax 
the criterion E has not yet reached the lower threshold 
e. Unfortunately, the conditions under which the algo- 
rithm terminates and returns a correct realizations are not 
known yet. 

Remark 3. (Convergence of Algorithm 3) . The criterion E 
of the algorithm always decreases or stays the same during 
the iterations, and the same holds for the cost functions in 
step 3-4. Hence, the criterion E and the cost functions in 
step 3-4 will converge to a fixed value. Whether this value 
is a global optimum remains a topic of future research. 
Note that the true parameters of a realization of / render 
E zero, but there might be many of them, see Example 1. 

5.1 Numerical example 

We conducted several numerical experiments. They were 
all performed on a PC with Intel Pentium 4 CPU 2.8 
GHz processor and 1.21 GB RAM memory under Windows 
XP, using Matlab 7.5 including the optimization toolbox. 
In line with the definitions above, we use K and N to 
quantify the dimension of the PWL of interest and M 
to denote the number of data points. Dynamical data 
{.f{l'i)jf{ti)}iLi was sampled from a given dynamical 
model X = f{x) at regular time intervals ti = i.A, i = 
1, . . . , M for some A > 0. Several initial states were chosen 
randomly, and a sampled trajectory was obtained for each 
of them. We then concatenated these sampled trajectories 
into one time series. Notice that concatenation of finite 
components of two output trajectories of a PWL is itself 
a finite component of a valid output trajectory. Hence, 
the combined data can be viewed as originating from the 
measurements of one valid output trajectory of a PWL. 
Gaussian noise with zero mean and variance was added 
to the obtained data. Finally, Algorithm 3 was applied to 
the resulting data. In the numerical experiments we first 
tested the performance of the algorithm on artificial PWL 
systems. Furthermore, we studied the algorithm on Lorenz 
systems and on the Tyson-Novak model for the cell cycle 
of budding yeast. 

Simulations on artificial PWL systems We generate an 
artificial PWL system containing K affine subsystems, 
each of the same dimension Uq ~ N. The state switching 
is obtained by partitioning the state space in K subsets 
{Vi, . . . ,Vk}- The switching is generated by using the 
following switching law: the discrete mode i is active, if 
the continuous state belongs to Vi. 

In the absence of noise, and if the subsystems are suffi- 
ciently sampled, the estimation algorithm is usually able 



Algorithm 3 Identification algorithm 
Inputs: data{/(t,),/(i,)}f£i 
Output: PWL S 



Fig. 1. 



Initialize the estimates of Aq,aq, q & Q, k :~ 0. 
repeat 

Solve the optimization for the weights 

w = arg mm^,^^^, J ^^^^^^ ^,^E{{ Aq,ag}qeQ:W* ) 

Vm = 1, . . . , M : ^ w;,„ = 1, V<z e g : wl„, G [0, 1] 
qeQ 

where for w = {wg.i}geQ,i=i....,Af, Aq e M"^", 
a, e M", g e g, 

1 

Note that the optimization problem above is a linear 
programming problem, hence the optimal values of 
w = {'Wq^i}q^Q^i^i^,,,^M take values in the set {0, 1}. 
Using the weights w = {w5,i}ggQ,i=i,...,Af recalcu- 
late the optimal values of Aq,aq, q G Q by solving 
the following minimization problem. 

{Aq,aq}qf,Q = arg min^, ^,_q^QE{{A*,a*q}qeQ,w) 

The optimization problem (5) is equivalent to the 
following linear least squares problem. Collect the 
unknowns Aq, aq, q £ Q into the matrix S 

S = [^1 ai A2 a2 ■ ■■ Ad an] 
Define the matrices L e Ri^+DDxAi ^^^j y e ]R"><*^ 

Y=[f{h) ■■■ f{tm)] 

wiiri ■ ■ ■ wmiTm 



_wiDri ■ ■ ■ wmotm 
n=[f{Uf l]^,Vi = l,...,M 

Then the solution of (5) can be read off from the 
entries of S, where 

S = arg min^g|j,„x(,.+i)D||i^ - SL\\^ 

5: Update E = E{{Aq,aq}qi^Q,w) using the new val- 
ues. 
6: until E < e 

7: Return S of the form (1) with Uq = n, Cq = /„, Cq = 
and Aq,aq as computed in the last iteration, for all 
qeQ, and Xq^o = M" for all q G Q. 



to find the system parameters, though occasionally it gets 
stuck in a local minimum. The amount of empirical data 
required to reconstruct the system parameters, within the 
bounds set by the noise, increases proportionally with N 
and K. Figure 1 below shows the result of the algorithm 
for N = 2, K ^ 5, M = 1800 data points. The Voronoi 
partitioning of state space of the original dataset (black) 
and the reconstructed dataset (red) are identical. The 
reconstructed system parameters and the original system 
parameters are identical as well. 




Fig. 2. 



Fig. 3. 





Further tests reveal that Algorithm 3 is able to find the 
correct parameters even if the regions {Vi, . . . ,Vk} are 
disconnected. 

When applying zero- mean Gaussian noise Af{0,a'^), the 
accuracy of the reconstruction decreases as a increases. 
Figure 2 below shows the result on a dataset of dimension 
A = 2 with K = 5 subsystems and M = 1800 data 
points with a SNR (signal-to-noise ratio) of 5%. Two 
subsystems could not be reconstructed, but the correlation 
between the three reconstructed systems with the best 
fitting original subsystems is 98.4 

Chaotic systems We applied Algorithm 3 to Lorenz 
systems Lorenz (1963). This example provides insight how 
the identification algorithm performs on systems that are 
not piecewise affine by nature. 

The Lorenz attractor can be visualized approximately as 
consisting of two linear subsystems Left and Right, with 
two different central attractors. A third region can be 
defined that involves the transition between the Left and 
Right subsystem. In Figure 3 the application of Algorithm 
3 to Lorenz systems with M = 5000 data points is 
presented for K = 2, 3. For K — 2 the two main parts 
of the Lorenz attractor are found. For K — 3 also the 
transition region from the left to the right attractor is 
found. Applying more than 3 subsystems does not improve 
the performance. 



Fig. 4. 
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Biological cell cycle models We applied Algorithm 3 to 
data simulated using the Tyson-Novak model Chen et al. 
(2004) for the yeast cell cycle. The cell cycle is the 
succession of events whereby a cell grows and divides 
into two daughter cells. The simulated data was used 
to construct a PWL approximation of the system using 
Algorithm 3. This resulted in a clear partition of the cell 
cycle dynamics as a function of the selected number of 
classes K. The best qualitative results were obtained for 
K = 2, see Figure 4. Moreover, for > 3 the results did 
not improve. It was found that the obtained two discrete 
modes correspond to the phases S (DNA synthesis) and 
M (mitosis) of the cell cycle. The reconstruction remains 
valid under noise with SNR at most 10%. 

6. CONCLUSIONS 

We presented some basic results on realization theory of 
PWLs and a practical identification algorithm. Analysis 
of the correctness and convergence of the presented algo- 
rithm remains a topic of future research. 
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